Stability analysis of sonic horizons in Bose-Einstein condensates 
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We examine the linear stability of various configurations in Bose-Einstein condensates with sonic 
horizons. These configurations are chosen in analogy with gravitational systems with a black hole 
horizon, a white hole horizon and a combination of both. We discuss the role of different boundary 
conditions in this stability analysis, paying special attention to their meaning in gravitational terms. 
We highlight that the stability of a given configuration, not only depends on its specific geometry, 
but especially on these boundary conditions. Under boundary conditions directly extrapolated from 
those in standard General Relativity, black hole configurations, white hole configurations and the 
combination of both into a black hole-white hole configuration are shown to be stable. However, 
we show that under other (less stringent) boundary conditions, configurations with a single black 
hole horizon remain stable, whereas white hole and black hole-white hole configurations develop 
instabilities associated to the presence of the sonic horizons. 

PACS numbers: 04.80.-y, 04.70.Dy, 03.75.Kk 
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I. INTRODUCTION 

It is widely expected that underneath the general rel- 
ativity description of gravitational phenomena there is 
a deeper layer in which quantum physics plays an im- 
portant role. However, at this stage we don't have 
enough intertwined theoretical and observational knowl- 
edge to know how an appropriate description of what 
underlies gravity is or should be. Moreover, starting 
from structurally-complete quantum theories of gravity, 
it could still be very difficult to extract the specific way in 
which the first "quantum" modifications to classical gen- 
eral relativity might show up (this happens for example 
within the Loop Quantum Gravity approach 

Analogue models of General Relativity (GR) 0,0] pro- 
vide specific and clear examples in which effective space- 
time structures ultimately emerge from (non-relativistic) 
quantum many-body systems. For certain (semiclassi- 
cal) configurations and low levels of resolution, one can 
appropriately describe the physical behaviour of the sys- 
tem by means of a classical (or quantum) field theory in 
a curved (Lorentzian) background geometry. However, 
when one probes the system with higher and higher res- 
olution, the geometrical structure progressively dissolves 
into a purely quantum regime Therefore, although 
analogue models cannot be considered at this stage com- 
plete models of quantum gravity (they do not lead to 
the Einstein equations in any regime or approximation), 
they provide specific and tractable models that repro- 
duce many aspects of the overall scenario expected in 
the realm of real gravity. 

The main objective of this and similar studies is to 
obtain specific indications about the type of deviations 
from the GR behaviour to be expected when quantum 
gravitational effects become important. All, under the 
assumption that the underlying structure to GR is some- 
what similar to that in condensed matter systems. In 



particular, in this paper we are interested in the be- 
haviour of gravity-like configurations containing horizons 
within Bose-Einstein condensates (BECs). (See, e.g., 
Refs. null] and H 1113 for reviews on BECs and for 
their usefulness as analogue models respectively). A nice 
feature of these systems is that their theoretical descrip- 
tion in terms of the Gross-Pitaevskii (GP) equation can 
be interpreted as incorporating, from the start, the first 
"quantum" corrections to the behaviour of the system. 
Linear perturbations over a background BEG configura- 
tion satisfy an equation which is a standard wave equa- 
tion over a curved effective spacetime plus corrections 
containing h. These corrections cause the dispersion re- 
lations in BEG to be "superluminal" (strictly speaking, 
supersonic): some perturbations can travel faster than 
the speed of sound in the system. The effects of these 
corrections in the linearized dynamical evolution of a 
configuration are especially relevant in the presence of 
horizons as their one-way-membrane nature simply dis- 
appears. This is in tune with the idea that a horizon can 
serve as a magnifying glass of the physics at high energies 
(see, e.g., Ref. [lOf). 

The specific objective of this paper is to analyze the dy- 
namical behaviour of (effectively) one-dimensional BEGs 
with density and velocity profiles containing one or two 
sonic horizons. In particular, we search for the presence 
of dynamical instabilities and analyze how their existence 
is related to the occurrence of these horizons. The stabil- 
ity analysis presented in Ref. foj for black hole- like config- 
urations with fluid sinks in their interior, concluded that 
these configurations were intrinsically unstable. How- 
ever, the WKB analysis of the stability of horizons in 
Ref. jll| suggested that black hole horizons might well 
be stable, while configurations with white hole horizons 
seem to posses unstable modes. Regarding configurations 
in which a black hole horizon is connected with a white 
hole horizon in a straight line, the analysis presented in 
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Ref. concluded, also within a WKB approximation, 
that these configurations were intrinsically unstable, pro- 
ducing a so-called "black hole laser" . However, when the 
white hole horizon is connected back to the black hole 
horizon to produce a ring, it was found that these con- 
figurations can be stable or unstable depending on their 
specific form |^ . This suggests that periodic boundary 
conditions can eliminate some of the instabilities associ- 
ated with the black hole laser. 

In this paper we will try to shed some light on all 
of these issues and clear up some of the apparent con- 
tradictions. To simplify matters, we will consider one- 
dimensional profiles that are piecewise uniform with ei- 
ther one or two step-like discontinuities. As has been 
discussed in Ref. ||9(], in terms of dynamical (in)stability, 
there seems to be no crucial qualitative difference be- 
tween the present case and a profile with smooth transi- 
tions between regions with an (asymptotically) uniform 
density distribution. Therefore the idealized case that we 
consider here should contain all the essential information 
relevant to more complicated profiles as well. The spe- 
cific way to examine the kind of instability we are inter- 
ested in, consists basically in seeking whether, under ap- 
propriate boundary conditions, there are complex eigen- 
frequencies of the system which lead to an exponential 
increase with time of the associated perturbations, i.e. 
a dynamical instability. Throughout the paper we will 
use a language and notation as close to GR as possible. 
In particular, we will use boundary conditions similar to 
those imposed in the standard quasi-normal mode anal- 
ysis of black holes in GR [T3| ■ One of the main results of 
our analysis is to highlight the fundamental importance 
of the boundary conditions in determining whether a con- 
figuration is stable or unstable. 

The structure of this paper is as follows. In the next 
section we will review the basic ingredients of gravita- 
tional analogies in BECs. At the same time, we will set 
up the conceptual framework for our discussion, based on 
a parametrization well adapted for an acoustic interpre- 
tation (a brief comparison with the Bogoliubov represen- 
tation is presented in appendix IB|| . Section IlIII contains 
a detailed formulation of our specific problem. This in- 
cludes the mode expansion in uniform sections, a deriva- 
tion of the matching conditions at each discontinuity and 
a discussion of the various boundary conditions to be ap- 
plied. Then, in section HVI we proceed case by case, an- 
alyzing different situations and presenting the results we 
have obtained for each of them. This includes a brief 
description of the numerical algorithm we have used. Fi- 
nally, in section [3 we discuss our results, compare them 
with other results available in the literature, and draw 
some conclusions. 



II. PRELIMINARIES 



yext(x)+g^'t*U, (2.1) 



ing the equation 



dt \ 2m 



where m is the boson mass, T4xt the external potential 
and g the coupling constant which is related to the cor- 
responding scattering length a through g — Airh^a/m. 
In this manner all quantum effects can, in principle, be 
taken into account. Once the Bose-Einstein condensation 
has taken place, the quantum field can be separated into 
a macroscopic wave function -0 (the corresponding order 
parameter) and a field operator (p describing quantum 
fluctuations over it: = ip + (p. The macroscopic wave 
function satisfies the Gross-Pitaevskii (GP) equation 



., 9 , , , 
dt 



2m 



v2-Hyext(x)-Hg|i^n V(i,x), 



(2.2) 



while for the linear quantum perturbation we have the 
Bogoliubov equation 

d f \ 
ih —p^ ( + ^cxt(x) -f g2|0n t^ + gV-' 

(2.3) 

Adopting the Madelung representation for the order pa- 
rameter 



w/h -\iit/h 



(2.4) 



(here n is the condensate density, the chemical poten- 
tial and a phase factor which is related to the velocity 
potential), and substituting in 1)2.2(1 we arrive at 

dtn = - — \7 ■ {nWB), (2.5a) 
m 

dtO = -7^(V0)2 -gn-V,^t-^^- Kjuantum, (2.5b) 

2m 

where the so-called "quantum potential" is defined as 

ft2 v^v^ 



quantum 



2m \/ri 



(2.6) 



In most situations the quantum potential in Eq. I|2.5bp 
can be neglected (see below). The resulting equations 



(|2.5a|l and l|2.5b|l are then equivalent to the continuity 
equation and the Bernoulli equation for a classical fluid. 
In this case, it is well known that the propagation of 
acoustic waves in the system can be described by means 
of an effective metric, thus providing the analo gy w ith the 
propagation of flelds in curved spacetimes [l^ll5l |. Given 
a background conflguration (no and Oq), this metric can 
be written as 



In second quantization, a dilute gas of interacting 
bosons can be described by a quantum field \l/ satisfy- 



-V 

1 



(2.7) 
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where (? = griQ/m. and v = WOo/m. These magnitudes, 
c and V, represent the local velocity of sound and the 
local velocity of the fluid flow respectively. 

The functions c{t, x) and v(i, x) completely charac- 
terize the acoustic metric. In GR any metric has to 
be obtained by solving the Einstein equations. Here, 
however, the magnitudes c(t, x) and v(t,x), and so the 
acoustic metric, are those satisfying the continuity and 
Bernoulli equations of hydrodynamics [Eqs. (|2.5() with- 
out the quantum potential]. Thus, these equations play 
a role analogous to the vacuum Einstein equations in 
GR. Of course, at the global non-linear level these equa- 
tions are completely different from the real Einstein equa- 
tions. But their way of acting when linearized around a 
background solution captures the essence of a proper lin- 
earized GR behaviour. 

There exist, however, situations in which the quantum 
potential in Eq. H2.5|l cannot be neglected. This is evi- 
dently the case if the characteristic length of the spatial 
variations of the condensate density is much smaller than 
the so-called healing length: f = H/{mc). But this case is 
not the only one. To illustrate this point, let us consider 
the dispersion law obtained for a homogeneous BEG (see 
below): 



the equations (|2.5(l to play the role of some sort of semi- 
classical vacuum Einstein equations. Their treatment is 
classical, but they incorporate corrections containing h. 
Therefore, BECs' standard treatment based on the GP 
equation provides an example of a way of incorporat- 
ing quantum corrections to the dynamics of a system 
without recurring to the standard procedures of back- 
reaction. Again, although at the global non-linear level 
these equations bear no relationship whatsoever with any 
sort of "semiclassical" Einstein equations, at the linear 
level, that is, in terms of linear tendencies of depart- 
ing from a given configuration, equations H2.5() encode 
the essence of the linearized GR behaviour (a Lorentzian 
wave equation in a curved background), semiclassically 
modified to incorporate a superluminal dispersion term, 
as we have already discussed. 

We will now proceed to describe the details of our spe- 
cific calculations. 



III. DYNAMICAL ANALYSIS 

As a first step in our calculations, let us linearize the 
Eqs. (|2.5|l . Let us write 



,2 7„2 



c'k 



(2.8) 



n(x, i) = no(x) +g ^ni(x, i), 
0(x,i) =^o(x)-f 0i(x,t), 



(3.1a) 
(3.1b) 



This is a "superluminally modified" dispersion relation 
due to the presence of the term with k*. For 1 we 

can rewrite this expression as 



where ni and 9i are small perturbations of the density 
and phase of the BEG. The Eqs. (|2.5(l then separate into 
two time-independent equations for the background. 



w ± c\/ 1 + ^C^/fc^ j k ~ {v±c)k+^cfk^+Oi^^k'^). 

(2.9) 



= -V • (c^v), 
1 



(3.2a) 



^2 y2^ 

= -^mv2-mc2-K=,t-M + ^ — , (3.2b) 



Here we clearly see that the (relative) importance of the pj^g time-dependent equations for the perturbations. 



term cx fc , given by 



£, k , depends not only on 



8(i>±c)'- 

the ratio between the corresponding wavelength and the 
healing length, ^k, but also on the specific features of the 
background magnitudes [note that the factor c/{v ± c) 
may be quite large]. In other words, the smallness of 
the corresponding wavelength is a necessary condition in 
order to neglect the contribution of the quantum poten- 
tial, but it is not a sufficient condition. One has to bear 
this issue in mind, especially when the system possesses 
horizons (i.e. points at which — v'^). 

Summarizing, there are background configurations 
which, when probed with sufficiently large wavelengths, 
act as if they were effective Lorentzian geometries. But 
there are other configurations for which this geometri- 
cal interpretation fails, irrespective of the probing wave- 
length. The latter situation occurs when there are hori- 
zons in the configuration: strictly speaking, we cannot 
talk about an "effective Lorentzian geometry" in the re- 
gions surrounding these horizons. 

Without forgetting this subtlety, we will continue to 
call (I2.7|l the "effective metric" in the system, even when 
analyzing the full GP equation. Then, we can consider 



dtui 
dtOi 



V • (niv 



^V0i) 



1 



(3.3a) 
(3.3b) 



We will restrict ourselves to work in (1-1-1) dimensions. 
This means that we consider perturbations propagating 
in a condensate in such a way that the transverse degrees 
of freedom are effectively frozen. In other words, the 
only allowed motions of both the perturbations and the 
condensate itself are along the x-axis. 

We will examine two types of one-dimensional back- 
ground profiles. The first type consists of two regions 
each with a uniform density and velocity, connected 
through a step-like discontinuity. The second type of pro- 
files consists of three homogeneous regions, and hence two 
discontinuities. We wish to know whether these profiles 
do or do not present dynamical instabilities. The un- 
derlying question is the relation between the presence of 
horizons and these dynamical instabilities. At each dis- 
continuity, matching conditions apply that connect the 
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magnitudes describing the condensate at both sides of 
the discontinuity. Furthermore, we need a set of bound- 
ary conditions, which determine what happens at the far 
ends of the condensate. Finally, in the uniform sections of 
the condensate, the regime can either be subsonic or su- 
personic, and of course there will be an acoustic horizon 
at each transition between a subsonic and a supersonic 
region. All these elements determine the characteristics 
of the system, and hence its eigenfrequencies. 

We will now describe these elements one by one in 
detail. 



A. Plane-wave expansion in uniform regions 

In order to study the dynamics of the system, let us 
first consider a region in which the condensate is homo- 
geneous (with c and v constant), and seek for solutions 
of Eqs. 1)3. 3|) in the form of plane waves: 

ni{x,t) = Ae'''^^-'^^\ (3.4a) 
ei{x,t) = Be''^'"=-'^'>\ (3.4b) 

where A and B are constant amplitude factors. Our aim 
is to elucidate about the possible instabilities of the sys- 
tem, so the frequency uj and the wavevector k in these 
expressions will be considered as complex hereafter [the 
existence of solutions with Im(a;) > would indicate the 
instability of the system]. Substituting into Eqs. H3.3|l 
we find 

/i{oj-vk) c^k'^ \ /A\ 

U 0. (3.5) 

For a non-trivial solution to exist, the determinant of the 
above matrix must vanish. This condition gives the dis- 



persion law (|2.8I) and since this is a fourth order polyno- 
mial in fc, its roots will, in general, give four independent 
solutions for the equations of motion in the form (|3.4|) . 



B. Matching conditions at a discontinuity 

Let us take a; = to be a point of discontinuity. The 
values of v and c both undergo a finite jump when cross- 
ing this point. These jumps have to satisfy the back- 
ground constraint vc^ = const [see Eq. (|3.2a|) ]. The so- 
lutions of Eqs. H3.3|l in the regions a; < and a; > have 
the form of plane waves which are then subject to match- 
ing conditions at a; = 0. It is not difhcult to see that 9i 
has to be continuous at the jump but with a discontin- 
uous derivative, while the function ni has to undergo 
a finite jump. The exact conditions can be obtained by 
integrating Eqs. (|3.3|l about an infinitesimal interval con- 
taining the point a; = 0. This results in the following four 
independent, generally valid matching conditions: 

[9-l] = 0, [vni + c^d^Oi] = 0, (3.6a) 




The square brackets in these expressions denote, for in- 
stance, [9i] = 6'i|^^Q+ — 6'i|^^g_. We can simphfy the sec- 
ond condition in Ea. ()3.6a|l to [c^(?2:6'i] = by noting that 
[vni] = because of the background continuity equation 
H3.2a|l , while for our choice of a homogeneous background 
the last condition becomes simply [9a;ni] = 0. 

For a given frequency tj, the general solution of Eqs. 
H3.3|l can be written as 



f 4 



"1 = ^8 



3=5 



icikj 



H^r—^e^''^-'^ (x>0). 



ic^kj 



(3.7) 



where {kj} are the roots of the corresponding dispersion equations (four roots for each homogeneous region), and 
the constants Aj have to be such that the matching conditions H3.6|) are satisfied. The subscripts L and R indicate 
the values of c and v in the left-hand-side (Ihs) and the right-hand-side (rhs) region respectively. We can write down 
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these conditions in matrix form AijAj = 0, where 



c?fc? 



1 

72- 



V fci 



1 

72- 



„2 1.2 

A,t- 



r -' k -' 



k4 
1 



kr, 
1 



r 



fee 
1 



1 



(3.8) 



Furthermore, these conditions have to be complemented 
with conditions at the boundaries of the system and then 
we will obtain the solution of a particular problem. 



C. Boundary conditions 

In order to extract the possible intrinsic instabilities of 
a BEC configuration, we have to analyze whether there 
are linear mode solutions with positive Im(w) that satisfy 
outgoing boundary conditions. By "outgoing" boundary 
conditions we mean that the group velocity is directed 
outwards (toward the boundaries of the system). The 
group velocity for a particular fc-mode is defined as 



Re 



duj 



Re 



c^k + lec^k^ 

LU — vk 



(3.9) 



where we have used the dispersion relation (|2.8() . The 
physical idea behind this outgoing boundary condition is 
that only disrupting disturbances originated inside the 
system can be called instabilities. 

To illustrate this assertion, let us look at the classi- 
cal linear stability analysis of a Schwarzschild black hole 
in GR. When considering outgoing boundary conditions 
both at the horizon and in the asymptotic region at in- 
finity p^ . only negative Im(a;) modes (the quasi-normal 
modes) are found, and thus the black hole configuration 
is stable. If the presence of ingoing waves at infinity were 
allowed, there would also exist positive Im(a;) solutions. 
In other words, the Schwarzschild solution in GR is sta- 
ble when considering only internal rearrangements of the 
configuration. If instead the black hole were allowed to 
absorb more and more energy coming from infinity, its 
configuration would continuously change and appear to 
be unstable. 

The introduction of modified dispersion relations adds 
an important difference with respect to the traditional 
boundary conditions used in linearized stability analy- 
sis in GR. Consider for example a black hole configura- 
tion. In a BEC black hole, one boundary is the stan- 
dard asymptotic region, just like in GR. For an acoustic 
(quadratic) dispersion relation, nothing can escape from 
the interior of a sonic black hole (the acoustic behaviour 
is analogous to linearized GR). But due to the super lu- 
minal corrections, information from the interior of the 



acoustic black hole can escape through the horizon and 
affect its exterior. Therefore, since we are taking this 
permeability of the horizon into consideration, the other 
boundary is not the black hole horizon itself (usually de- 
scribed in GR by an infinite value of the "tortoise" coor- 
dinate; see for example 1^), but the internal singularity. 
The outgoing boundary condition at such a singularity 
reflects the fact that no information can escape from it. 

There is another complication that deserves some at- 
tention. In the case of a phononic dispersion relation, the 
signs of Wg = V ± c and Im(fc) = (v ± c)Im(ti;) coincide 
for Im(aj) > 0. For example, in an asymptotic x — > -foo 
subsonic region, an outgoing fc-mode has Vg > 0, so that 
Im(fc) > and, therefore, the mode is damped towards 
this infinity (giving a finite contribution to its norm). 
Owing to this fact, in the linear stability analysis of 
black hole configurations, it is usual to assume that sta- 
ble modes correspond to non-normalizable perturbations 
(think of the standard quasi- normal modes), while un- 
stable modes correspond to normalizable perturbations. 

When considering modified dispersion relations, how- 
ever, this association no longer holds. In particular, 
with the BEC dispersion relation, in an asymptotic 
X — > -|-oo region, among the unstable (lm{uj) > 0) out- 
going (vg > 0) fc-modes, there are modes with Im(fc) > 
as well as modes with Im(fc) < 0. An appropriate in- 
terpretation of these two possibilities seems to be the 
following. The unstable outgoing modes that are con- 
vergent at infinity (those with Im(fc) > 0) are associated 
with perturbations of the system that are initiated in an 
internal compact region of the system. Unstable outgo- 
ing modes that are divergent at infinity are associated 
with initial perturbations acting also at the boundary at 
infinity itself. 

Take for example a black hole-like configuration of the 
form described in Fig. ^ The right asymptotic region, 
which can be interpreted as containing a "source" of 
BEC gas in our analogue model, simulates the asymp- 
totic infinity outside the black hole in GR. The conver- 
gence condition at the rhs then implies that the perturba- 
tions are not allowed to affect this asymptotic infinity ini- 
tially. However, for the left asymptotic region, this con- 
dition is less obvious. In our BEC configuration, this left 
asymptotic region can be seen as representing a "sink" . 
It corresponds to the GR singularity of a gravitational 
black hole. The fact that in GR this singularity is sit- 
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FIG. 1: Flow and sound velocity profiles with step-like dis- 
continuities simulating a black hole-like configuration. The 
negative value of v indicates that the fluid is left-moving. At 
the rhs, the fluid is subsonic since c > |?;|. At the Ihs it has 
become supersonic. At a; = 0, there is a sonic horizon. 

uated at a finite distance (strictly speaking, at a finite 
amount of proper time) from the horizon, indicates that 
it might be sensible to allow the perturbations to affect 
this left asymptotic region from the start. We will there- 
fore consider two possibilities for the boundary condition 
at a; ^ — oo. (a) Either we impose convergence in both 
asymptotic regions, thereby eliminating the possibility 
that perturbations have an immediate initial effect on the 
sink, or (b) we allow the perturbations to affect the sink 
right from the start, i.e. we don't impose convergence 
at the left asymptotic region. The option of imposing 
the convergence at the left asymptotic region could be 
interpreted as excluding the influence of the singularity 
on the stability of the system. In other words, condition 
(b) would then be equivalent to examining the stability 
due to the combined influence of the horizon and the sin- 
gularity, while under condition (a) only the stability of 
the horizon would be taken into account. 

As a final note to this discussion, since we are inter- 
ested in the analogy with gravity, we have assumed an 
infinite system at the rhs. In a realistic condensate other 
boundary conditions could apply, for example taking into 
account the reflection at the ends of the condensate (see 
e.g. Hi). 

IV. CASE BY CASE ANALYSIS AND RESULTS 

We will now briefly describe the general calculation 
method which wc have used, and then discuss case by 
case the specific configurations we have analyzed. 

A. Numerical method 

We first consider background flows and sound velocity 
profiles with one discontinuity. We will always assume 
left-moving flows. 

We are seeking for possible solutions of the linearized 
Eqs. with Im((jj) > 0. We use the following numer- 

ical method: 

1. For each frequency w in a grid covering an appro- 



priate region of the upper-half complex plane, we 
calculate its associated fc-roots [by solving the dis- 
persion relation (|2.8|l ] and their respective group 
velocities (|3.9|) at both the Ihs and the rhs of the 
configuration. 

2. We then take the four equations AijAj — 0, where 
A is the 4x8 matrix H3.8|l determined by the match- 
ing conditions at the discontinuity. For each mode 
kj that does not satisfy the boundary conditions 
in the relevant asymptotic region, we add an equa- 
tion of the form Aj =0. Thus we have a total set 
of equations which can be written as ^ijAj = 0, 
where A is now a (4 -I- N) x 8 matrix the number N 
of forbidden modes can in principle vary between 
and 8). Numerically it is convenient to normalize 
A in such a way that its rows are unit vectors. 

3. We can then define a non-negative function F{uj), 
where F[u;) = means that the frequency to is an 
eigenfrequency of the system, in the following way. 

• If iV < 4, then F{uj) = 0. Indeed, we have 8 
variables Aj and 4 -I- TV equations. Then, it 
is obvious that there will always exist a non- 
trivial solution {Aj}. 

• If iV = 4, then F{uj) = |det(A)|. In this case 
there will be a non-trivial solution only if the 
determinant of the matrix A vanishes. 

• If = 5, then F{oj) is taken to be the sum of 
the modulus of all possible determinants that 
are obtained from A by eliminating one row. 
Notice that in this case A is a non-square 9x8 
matrix because there are more conditions than 
variables. In this situation it is highly unex- 
pectable to find zeros in F as this would mean 
a double degeneracy. 

• If > 5, F{uj) is defined by a straightforward 
generalization of the procedure for N — 5. 

4. We plot the function F{uj) in the upper half of the 
complex plane, and look for its zeros. Each of these 
zeros indicates an unstable eigenfrequency, and so 
the presence (or absence) of these zeros will indicate 
the instability (or stability) of the system. 

In all our numerical calculations we have chosen val- 
ues for the speed of sound and the fluid velocity close 
to unity. Moreover, we set vc^ — 1 and choose units such 
that ^ c = 1. The typical values of the velocity of sound in 
BECs range between lmm/s~10mm/s, while the healing 
length lies between 10~^mm - 10~^mm. In consequence, 
our numerical results can be translated to realistic phys- 
ical numbers by using nanometres and microseconds as 
natural units. For example, the typical lifetime for the 
development of an instability with Im([jj) ~ 0.1 would be 
about 10 microseconds. We have checked that our re- 
sults do not depend on the particular values chosen for 
the velocities of the system. 
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B. Black hole configurations 

Consider a flow accelerating from a subsonic regime 
on the rhs to a supersonic regime on the Ihs, see Fig- 
ure n For rhs observers this configuration possesses a 
black hole horizon. For such configurations with a sin- 
gle black hole-like horizon, when requiring convergence 
in both asymptotic regions [case (a)], there are no ze- 
ros (see Figure EJ, except for two isolated points on the 
imaginary axis (see Figure |31 which is a zoom of the rel- 
evant area in Figure [2Jl- (We always check the existence 
of a zero by zooming in on the area around its location 
up to the numerical resolution of our program.) These 
points are of a very special nature. They are located at 
the boundary between regions with different number N 
of forbidden modes in the asymptotic regions. The zeros 
that we will find for other configurations are of a totally 
different nature: they are sharp vanishing local minima 
of F{uj) living well inside an area with a constant value 
of N (N = 4 to be precise). We discuss the meaning of 
these special points in Appendix^ For now, let it suf- 
fice to mention that points of this kind are always present 
in any flow, independently of whether it reaches super- 
sonic regimes or not. Hence it seems that they do not 
correspond to real physical instabilities, since otherwise 
any type of flow would appear to be unstable. Accord- 
ingly, in the following, we will not take these points into 
consideration. When we assert that a figure is devoid of 
instabilities, we will mean that the function F{ll}) has no 
zeros except for the special ones just mentioned. 

Figure |21 also shows that the system remains stable 
even when eliminating the condition of convergence at 
the Ihs [case (b)]. 

To sum up, configurations possessing a (single) black 
hole horizon are stable under the general boundary 
conditions that we have described, i.e. outgoing in 
both asymptotic regions and convergent in the upstream 
asymptotic region, independently of whether convergence 
is also fulfilled in the downstream asymptotic region or 
not. 



C. White hole configurations 

Let us now consider flows decelerating from a super- 
sonic regime (rhs) to a subsonic one (Ihs). From the point 
of view of Ihs observers, the geometric conflguration pos- 
sesses a white hole horizon. In GR, a white hole corre- 
sponds to the time reversal of a black hole. Therefore, 
unstable modes of a white hole configuration would cor- 
respond to stable (quasinormal) modes of the black hole. 
However, when modified dispersion relations are present, 
the precise definition of quasinormal modes cannot be 
based only upon the outgoing character of the modes, 
but their divergent or convergent character also has to be 
taken into account. Having in mind that in the acoustic 
approximation (proper Lorentzian behaviour) the outgo- 
ing character of a quasinormal mode implies that this 



mode diverges at the boundaries at infinity, it is rea- 
sonable to impose divergence as an additional defining 
requirement (apart from being outgoing) for a quasinor- 
mal mode in the presence of modified dispersion rela- 
tions. Using this definition, we have checked (by analyz- 
ing the lower- half complex co plane), that black hole con- 
figurations do not show any quasinormal (stable) eigen- 
frequency. Thus, we can conclude that one- dimensional 
white holes are stable. We emphasize here the word one- 
dimensional because we do not expect this situation to 
remain true in higher dimensions. We know for example 
that standard GR black holes in 3-1-1 dimensions pos- 
sess quasinormal modes. We expect these quasinormal 
modes to subsist when taking into account departures 
from the acoustic (Lorentzian) dispersion relation; we 
only expect them to acquire modified eigenfrequencies. 
These quasinormal modes would then identify instabili- 
ties of the corresponding white hole configuration. We 
leave the analysis of the quasinormal modes in different 
analogue gravitational configurations in BECs for future 
work, since this analysis has its own subtleties. 

The boundary conditions appropriate for the analysis 
of white hole- like configurations correspond to only hav- 
ing ingoing waves (due to time reversal) at the bound- 
aries. But from the point of view of acoustic models in 
a laboratory, the analysis of the intrinsic stability of the 
flow (under the outgoing boundary conditions described 
above) is also interesting. This analysis also has par- 
ticular relevance with regard to configurations with two 
horizons (see below). 

In Fig. [51 we see that under outgoing boundary condi- 
tions the flow is stable when convergence is required in 
both asymptotic regions [case (a)], but exhibits a con- 
tinuous region of instabilities at low frequencies when 
convergence is fulflUed only at the rhs [case (b)] . Indeed, 
in this continuous region = 3, in other words the al- 
gebraic system KijAij = is underdetermined and any 
frequency is automatically an eigenfrequency. 

When looking at the case of a completely subsonic flow 
suffering a deceleration (see Fig. we find something 
similar. The system is clearly stable when convergence is 
imposed at the Ihs, i.e. in case (a). Without convergence 
at the Ihs, case (b), there is a continuous strip of insta- 
bilities which corresponds, as in the white hole case, to 
a region where A^ = 3. However, this region is localized 
at relatively high frequencies and so disconnected from 
w ~ 0. We can say that part of the continuous region of 
instabilities found in the white hole configuration has its 
origin merely in the deceleration of the fiow (giving rise 
to this high frequency strip). But there is still a complete 
region of instabilities that is genuine of the existence of 
a white hole horizon. In fact, by decreasing the healing 
length parameter ^, the strip moves up to higher and 
higher frequencies, becoming less and less important as 
one approches the acoustic limit. However, the contin- 
uous region of instabilities associable with the horizon 
does not change its character in this process. 

We can therefore conclude the following with regard 



Black hole White hole 

Csupcr = 0.7, Csub = 1-8. Csub = 1.8, Csupcr = 0.7. 

JV(ai) log[-P(«i)] iog[F(m)] 




-4 -2 2 4 -4 -2 2 4 -4 -2 2 4 -4 -2 2 4 

Re(ra) Re(ffl) Re (to) Re (to) 

FIG. 2: Stability analysis under outgoing boundary conditions for profiles with one discontinuity. Represented is the relevant portion of the upper-half complex 
frequency plane. From top to bottom and left to right: black hole and white hole configurations, accelerating and decelerating subsonic flows (the speed of sound c 
is indicated for each region and the velocity v is then obtained from the constrain vc^ = 1; c > 1 corresponds to a subsonic region, c < 1 to a supersonic region; in 
addition, we use = 1 in all our calculations). The Ihs pictures represent the number TV of forbidden modes in the asymptotic regions. The rhs pictures represent 
the function F{u!) (to enhance the contrast, we have drawn the logarithm), and white points or regions, where F{u>) = 0, represent instabilities. In the upper pictures 
[case (a)], convergence has been imposed in both asymptotic regions. In the lower pictures [case (b)], convergence has been imposed only in the upstream asymptotic 
region. It is seen that black hole configurations are stable in both case (a) and (b), as are accelerating subsonic fiows. White hole configurations are stable in case (a), 
but develop a huge continuous region of instabilities in case (b). Only a small strip of instabilities subsists in the decelerating subsonic flow, indicating that the major 
part of this unstable region is a genuine consequence of the existence of the white hole horizon. Note that continuous regions of instability correspond to < 4. 
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FIG. 3: Two special zeros of the function F(u) appear in the 
stability analysis of a black hole configuration (this plot is a 
zoom of the corresponding plot in Fig. I^J. They are located 
at the boundary between regions with different number A'^ of 
prohibited modes. These points do not seem to represent real 
instabilities of the system (see appendix 0. 

to decelerating configurations. When convergence is ful- 
filled downstream, the configuration is stable, regardless 
of whether it contains a white hole horizon or not. When 
this convergence condition is dropped, there is a tendency 
to destabilization. In the presence of a white hole hori- 
zon, the configuration actually becomes dramatically un- 
stable, since there is a huge continuous region of instabil- 
ities, and even perturbations with arbitrarily small fre- 
quencies destabilize the configuration. In the absence of 
such a horizon, only a small high-frequency part of this 
unstable region subsists. 



D. Black hole— white hole configurations 

Consider flows passing from being subsonic to super- 
sonic and then back to subsonic (Figure The numer- 
ical algorithm we have followed to deal with this prob- 
lem is equivalent to the one presented above, but with a 
larger set of equations. In this case we have 12 arbitrary 
constants Aj, which have to satisfy 8 + equations: 4 
matching conditions at each discontinuity and A^(0 — 8) 
additional conditions of the form Aj — Q, corresponding 
to modes that do not fulfill the boundary conditions in a 
particular asymptotic region. 

When convergence is imposed at the Ihs, we do not 
find any instabilities, regardless of whether the fiuid is 
globally accelerating or decelerating [the final Ihs fluid 
velocity is larger or smaller than the initial rhs one re- 
spectively, see Figure |S1 cases (a)]. Also when replac- 
ing the intermediate supersonic region by a subsonic one, 
thereby removing the acoustic horizons, the fluid is sta- 
ble, independently of whether it is globally accelerating 
or decelerating. 



x=-L/2 






c 

x=L/2 ? 










V 


-c 



FIG. 4: Flow and sound velocity profiles with step-like dis- 
continuities simulating a black hole-white hole configuration. 



When dropping the convergence condition at the Ihs 
the situation changes completely. When the intermedi- 
ate region is supersonic, i.e. in a black hole-white hole 
configuration, a discrete set of instabilities appears at 
low frequencies [Fig. [peases (b)]. It is worth mentioning 
that, when carefully looking at plots of type (a)-cases, 
we observe some traces of these zeros in the form of lo- 
cal minima which can be understood as particularly soft 
regions. These regions, although very close to zero in 
some situations, never give rise to real zeros, as we have 
carefully checked by zooming in. Notice that these local 
minima appear in regions with TV = 5 where a zero would 
mean a double degeneracy within the row vectors in the 
corresponding matrix Ay . When the fiuid is globally de- 
celerating, additionally there is a continuous region of 
instabilities at higher frequencies. Indeed, in this region, 
as in the case of the white hole configuration, < 4, and 
so every frequency in this region automatically represents 
an instability. When the intermediate region is subsonic, 
the discrete set of local minima at low frequencies disap- 
pears, but the continuous strip of instabilities at higher 
frequencies persists in the case of a globally decelerat- 
ing fluid. The discrete set of instabilities is therefore a 
genuine consequence of the existence of horizons. 

The number of discrete zeros we find in the black hole- 
white hole configuration increases with the size L of the 
supersonic region (see Fig. EJ, while their Im(w) de- 
creases. This suggests that the region between the hori- 
zons acts as a sort of well discretizing some of the in- 
stabilities found for the white hole configurations. The 
larger the well, the larger the amount of instabilities, but 
the longer- lived these instabilities. 

To summarize, when requiring convergence in both 
asymptotic regions, all the types of configurations with 
two discontinuities that we have discussed are stable. 
When not requiring convergence at the Ihs, discretized 
instabilities appear associated with the presence of hori- 
zons. 



E. Black hole configurations with modified 
boundary conditions 

We have seen in section HV Bl that configurations with 
a single black hole horizon do not possess instabilities in 



(a) 



(b) 
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Black hole— white hole (globally accelerating) 

Csub-lhs = 1-8, Csupcr = 0.7, Cgub-rhs = 1-9. 
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-2 
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Globally accelerating subsonic flow 

Csub-lhs — l'^, Csub — 1.8, Csub-rhs — 1.9. 

l0g[i'(ffl)] 



N(as) 




n 
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Black hole— white hole (globally decelerating) 

Csub-lhs = 1.9, Csupcr ~ 0.7, Csub-rhs ~ 1.8. 

log[/-(ra)] 




□ 



-4-2 



Re(c>) 



Globally decelerating subsonic flow 



Csub-lhs — 1.9, Csub 
N(m) 



1.8, Csub-rhs ~ 1.7. 

l0g[i'(ffl)] 





Re (lis) 



Re(ffl) 



FIG. 5: Stability analysis for profiles with two discontinuities. From top to bottom and left to right: globally accelerating and decelerating black hole-white hole 
configurations, globally accelerating and decelerating subsonic configurations. In all these plots we have used L = 2.5 as the size of the intermediate region (see 
also caption under Fig. |2j. When convergence is imposed in both asymptotic regions [case (a)], all the configurations are stable. When convergence is only imposed 
upstream [case (b)], the configurations with sonic horizons present a discrete set of instabilities at low frequencies, while the decelerating configurations show a small 
continuous unstable strip at high frequencies. The decelerating configuration with sonic horizons combines both types of instabilities. 
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FIG. 6: The number of discrete instabilities in the black hole- 
white hole configuration increases with the size L of the in- 
termediate supersonic region. In these plots we have used 
Caub-ihs = 1.8, Cauper = 0.7, Csub-rhs — 1.9 (with their corre- 
sponding V = 1/c^) and ^ c = 1. 



any situation. However, as we have just discussed, when 
they are continued into a white hole configuration, some 
instabiUties can show up. Here, we would like to point 
out that the same happens if instead of extending the 
black hole configuration we introduce a wall (or sink) at 
a finite distance inside the supersonic region, described 
by other boundary conditions than the ones we have con- 
sidered so far. For example, by replacing the Ihs bound- 
ary conditions by 9\x=-l — 0, we obtain Fig. |7| We 
can perfectly see how a set of discrete unstable modes 
appears. 



V. DISCUSSION AND CONCLUSIONS 

Let us start by discussing the stability of configura- 
tions with a single black hole-like horizon in analogue 



log[i?(co)] 



-2-1012 

Re(cD) 

FIG. 7: A discrete set of instabilities appears in a black hole 
configuration when the Ihs asymptotic region representing the 
singularity is replaced by a wall or sink. In this plot we have 
used Csub = 1.9 and Csuper = 0.7 (with their corresponding 
V = 1/c^); in addition we have taken L = 6 as the size of the, 
now finite, internal region. 



systems that incorporate superluminal dispersion rela- 
tions. We have seen that by requiring purely outgoing 
and convergent boundary conditions in both asymptotic 
regions, these configurations do not show any signs of 
instability. The same applies when dropping the con- 
vergence condition downstream (i.e. on the Ihs). This 
seems to contradict the results in Ref. @ . There, the ex- 
istence of a future (spacelike) singularity inside the black 
hole, from which no information is allowed to escape, was 
implemented by introducing a sink in the supersonic re- 
gion at a finite distance from the horizon. Then, it was 
found that there were discrete instabilities in the system. 
However, these instabilities correspond to the following 
particular set of boundary conditions: i) At the asymp- 
totic region, only convergent boundary conditions were 
imposed, without any condition about the direction of 
propagation (in- or outgoing) of the perturbations; ii) At 
the sink, two types of boundary conditions were required, 
specifically designed for dealing with symmetric and anti- 
symmetric configurations. In our language these bound- 
ary conditions correspond to {0'|sink = 0, n'|sink = 0} and 
{^Isink = 0, n|sink = 0} respectively. In comparing this re- 
sult with ours we have checked two important facts. On 
the one hand, their unstable modes have ingoing contri- 
butions at the asymptotic region. On the other hand, 
the boundary conditions at the sink are such that they 
combine outgoing and ingoing contributions - their sink 
implementation makes waves reaching the sink bounce 
back towards the horizon. These two facts are respon- 
sible for the unstable behaviour of these black hole-like 
configurations. If no energy is introduced into the sys- 
tem from the asymptotic region (in other words, if only 
outgoing perturbations are allowed) and moreover any 
bouncing at the sink is eliminated, then these configu- 
rations are stable. This is in agreement with the result 
found in Ref. .11.] . 

In the case of configurations with a single white hole 
horizon, we have seen that with outgoing and conver- 
gent boundary conditions in both asymptotic regions, 
there are no instabilities in the system. However, when 
eliminating the convergence condition in the downstream 
asymptotic region, one finds a continuous region of insta- 
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bilities surrounding uj — 0. Thus, we see that these white 
hole configurations are stable only when the boundary 
conditions are suSiciently restrictive. 

When analyzing configurations connecting two difi^er- 
ent subsonic regions, we have also seen that, again, 
when convergence is required at the Ihs, they are stable. 
But when this convergence condition is relaxed, glob- 
ally decelerating configurations tend to become unstable, 
whereas globally accelerating ones remain stable. The 
instabilities of these decelerating configurations without 
horizons (i.e. purely subsonic ones) show up, however, in 
a small strip at high frequencies. In contrast, white hole 
configurations present instabilities for a wide range of 
frequencies, starting from arbitrarily small values. This 
points out that the presence of a white hole horizon dras- 
tically stimulates the instability of the configuration. 

With regard to the black hole-white hole configura- 
tions, we have seen that, as before, with outgoing and 
convergent boundary conditions, they are stable. How- 
ever, when relaxing the convergence condition down- 
stream, they develop a discrete set of unstable modes. 

In the analysis of the black hole laser instability in 
Ref. the authors found that these black hole-white 
hole configurations were intrinsically unstable. However, 
they did not analyze what happens to the modes at 
the Ihs infinity. Our analysis shows that by restricting 
the possible behaviour of the modes in the downstream 
asymptotic region, one can eliminate the unstable be- 
haviour of the black hole laser. This is in agreement 
with the results in Refs. There, the instabilities 

can in some cases be removed by requiring periodicity, 
that is, by imposing additional boundary conditions to 
the modes. 

To sum up, wc have shown the high sensibility of the 
stability not only on the type of configuration (the pres- 
ence of a single horizon or of two horizons, the acceler- 
ating or decelerating character of the fluid) , but particu- 
larly on the boundary conditions. With outgoing bound- 
ary conditions, when requiring convergence at the down- 
stream asymptotic region, both black hole and white hole 
configurations are stable (and also the combination of 
both into a black hole-white hole configuration). When 
relaxing this convergence condition at the Ihs, configu- 
rations with a single black hole horizon remain stable, 
whereas white hole and black hole-white hole configura- 
tions develop instabilities not present in (subsonic) flows 
without horizons. 
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APPENDIX A: ZEROS AT THE BOUNDARIES 
OF THE REGIONS IN N{lo) 

Given an w, one can flnd its four associated k roots, 
{kj}. If instead of co one takes (D = — w*, it can be seen 
that the new roots {kj} arc just {— fc*}. For this reason 
the function F(uj) is mirror symmetric with respect to 
the imaginary axis (this is seen in all our flgures). Now, 
when to is pure imaginary (to = —lo*), the set {kj} has to 
be equal to the set {—k*}. There are three posibilities. 
Either all four roots are pure imaginary, two are imag- 
inary and the other two complex satisfying kj = — fc^*, 
with j ^ I, or there are two pairs of complex roots sat- 
isfying kj = —k^. When moving through the imaginary 
UJ axis, there are points at which there is a transition 
from one of these possibilities to another. At any transi- 
tion point there has to be a pair of imaginary roots with 
equal value. Deflning to" = Im(w) and k = —ik, the 
dispersion equation H2.8|l can be written as 

_ y^f _ (^c^ „ jc^f = 0. (Al) 

This is a fourth order polynomial in k with real coeffi- 
cients. If this polynomial has two equal real roots then 
we know that the derivative with respect to k of the poly- 
nomial has to be zero at that point. It is not difficult to 
see that this also implies that the derivative with respect 
to K of the function 

/ 1 X 1/2 

{w -vk)t ic^ - ^c^^^K^j K (A2) 

has to be zero at this same point. But this derivative 
coincides with the definition of the group velocity given 
in l|3.9|l (when uj and k are pure imaginary, duj/dk is 
directly real). Therefore, we conclude that at any transi- 
tion point on the imaginary w-axis we have degeneracy: 
at least two imaginary k roots with equal value. At the 
same time, the group velocity associated with them be- 
comes zero. This is why these points are located at the 
boundary between regions with a different number of for- 
bidden modes: these are places in which outgoing modes 
transform into ingoing ones. The zero that appears in 
the function F{ll}) at these points is due to the degener- 
acy and does not tell us anything about the existence or 
not of a real instability there. To know whether a real 
instability appears, one first has to find the actual four 
independent solutions of equations (|3.3|l at the point that 
led to the degeneracy. Let us check under which circum- 
stances one can find a solution of the form 

ni(2;,t) = Aixe'C^^-'^)*, (A3a) 
Oi (x, t) = Bix e*(fe^~'^)* + ^^e^C^^""')* . (A3b) 

For these expressions to be a solution of Eqs. (|3.3|l . the 
following conditions have to be satisfied: 
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i{u} — vk) 



21,2 



c'k 



= 0, 



l + \(^k^ -i{uj-vk)) \Bi^ 



(A4) 



i[uj — vk) c^k"^ 
1 + -iiuJ-vk) 




-2ic^k\ I Ax 




0. 



(A5) 



From the first condition we obtain that the dispersion 
relation (|2.8|l has to be fulfiUcd. As a consequence we 
also find that -Bi = A\(^u) — vk) / [ic^k"^). Now, from the 
second condition we obtain 

^[u~vk)j^ \_i,^2j^ , ] 

(A6) 

This is a system of two equations from which, eliminating 
A1/B2 and after some rearranging, we obtain: 

c^k+]^c^^'^k'^ + v{oj-vk)^0. (A7) 

This is exactly the condition for a vanishing group ve- 
locity H3.9|l . Therefore, when functions in the form of 
plane waves do not lead to four linearly independent so- 
lutions, but for example two are "degenerate" , then we 
can use the previous solution HA3|I avoiding this degener- 
acy. Once we have the actual four independent solutions 
of the problem, they have to be matched with the four 
solutions in the other region (typically these will have 
the form of plane waves, unless we are in a very special 
situation in which degeneracy occurs in both regions at 
the same time) and see whether there is a combination 
satisfying all the boundary conditions. 

Although we haven't made such a full detailed calcu- 
lation, the fact that this kind of situation occurs in any 
type of flow indicates that it is safe to assume that they 
do not represent real instabilities, as already mentioned. 



APPENDIX B: BOGOLIUBOV 
REPRESENTATION 

All calculations and numerical simulations presented 
in the text have been performed independently by using 
the acoustic representation and the Bogoliubov represen- 
tation described in this appendix. We have found 
identical results with the two methods, double checking 
in this way the absence of numerical artifacts. 

Consider a one-dimensional setup with a potential V^xt 
that produces a profile for the speed of sound of the form 



{Co, a: < 

co[l + (o- - l)a;/e], 0<x<e. (Bl) 
CTCo, t < X 

We will assume a > \ and a flow velocity in the inward 
[x — > —00) direction, i.e. a black hole-like configuration 
for a rhs observer. The limit e ^ provides the same 
profile we have discussed in the main text. 

The condensate wave function can be written as the 
sum of a stationary background state '00 a-nd a pertur- 
bation (j) satisfying 

idt(t> = -dl<j)l2 -f (c^ - w 2/2 -f 9^c/2c)0 -f c^e^' , 

(B2) 

that will then be expanded into Bogoliubov modes 

ct>^u^{x)e-'^' +wl{x)e'^'K (B3) 

^From now on, we drop the subindex w, but remember 
that all equations should be valid for every (complex) 
frequency u separately. 

The assumption of a small e leads to a linear solu- 
tion for the modes u and w in the intermediate region 
< a; < e. Together with the transition conditions at 
a; = and x — t, which substitute the singular charac- 
ter of d^cjc at those points, this leads to the following 
connection formulas (in the limit e ^ 0), see Ref. 0: 

[w/c] =0, \d^{cu)\ =0, (B4) 

\wlc\ =0, \d^{cw)\ =0, (B5) 

where as before we have used the notation 

Let us write the density and phase in terms of the 
stationary values plus perturbations: 

n-no+g-iHi, = 00 + ^1, (B6) 

and note that the matching condition for the background 
velocity potential is — 0. 

The relation between the acoustic and the Bogoliubov 
representation of the perturbations can easily be ob- 
tained by noting that, to first order, 

^ .J^gc{n^/2mc^ + 161)^^" , (B7) 
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where we have used the fact that = gno/m. Note that 
ni and 9i are real. Then, comparison of Eq. l|B7p with 
the Bogohubov mode expansion = ite~*'^* + w*e*" * 
yields 

u = v/m/gc(ni/2mc^ + 16*1)6*^", (B8) 
w = ^/mjgc{ni/2mc^ ~ i9i)e-''''>. (B9) 

We then have 

[u/c] cx [{ni/2c^ + iei)e'''"] cx [ni/2c'^ + i0i] (BIO) 

and likewise for w. Comparing with Eq. (|B4|I . we obtain 

[ni/c^] cx [u/c] + [w/c] = 0, (Bll) 
[6*1] cx [u/c] - [w/c] ^ 0, (B12) 



which correspond to two of the conditions in Eqs. 

Taking into account that [9q] — [c^v] = [9i] = and 
that dxCL = dxCR = 0, we can write 

[dxicu)] oc i[vni]/2 + [dxni]/2 + i[c'dxei], (B13) 
[dx{cw)] cx - i[vni]/2 + [dxni]/2 ~ i[c'dxei]. (B14) 

In other words [9a;«^i] = and [vni/2 + c^dxOi] — 0. Be- 
cause of the continuity equation vc^ = const, we have 
[vni] cx [ni/c^] — 0, so that the boundary conditions ul- 
timately become 

[dxui] = 0, [c^OxOi] = 0, (B15) 

i.e. the other two conditions in Eqs. H3.6|) with the cor- 
responding simplifications. 
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